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Abstract 

A combination of analytic modeling and sys- 
tem identification methods have been used to de- 
velop an improved dynamic model describing the 
response of articulated rotor helicopters to control 
inputs. A high-order linearized model of coupled 
rotor/body dynamics including flap and lag degrees 
of freedom and inflow dynamics with literal coeffi- 
cients is compared to flight test data from single ro- 
tor helicopters in the near hover trim condition. 
The identification problem was formulated using 
the maximum likelihood function in the time do- 
main. The dynamic model with literal coefficients 
was used to generate the model states, and the mod- 
el was parametrized in terms of physical constants 
of the aircraft rather than the stability derivatives, 
resulting in a significant reduction in the number of 
quantities to be identified. The likelihood function 
was optimized using the genetic algorithm ap- 
proach. This method proved highly effective in 
producing an estimated model from flight test data 
which included coupled fuselage/rotor dynamics. 
Using this approach it has been shown that blade 
flexibility is a significant contributing factor to the 
discrepancies between theory and experiment 
shown in previous studies. Addition of flexible 
modes, properly incorporating the constraint due to 
the lag dampers, results in excellent agreement be- 
tween flight test and theory, especially in the high 
frequency range. 


Presented at Piloting Vertical Flight Aircraft: A 
Conference On Flying Qualities and Human Fac- 
tors, San Francisco, California, 1993. 


Introduction 

The investigation of rotorcraft dynamics, and 
specifically the coupled fuselage/rotor dynamics, is 
motivated by increasing sophistication in rotorcraft 
stability analyses and by the emergence of high- 
performance flight control system design require- 
ments. The past few years have seen a concentrated 
effort directed toward providing an analytic simula- 
tion model of coupled fuselage/rotor dynamics and 
model validation against flight test data. 

Helicopter dynamics include the rigid-body 
responses demonstrated by fixed-wing aircraft, 
plus higher-frequency modes generated by the in- 
teractions of the rotor system with the fuselage. For 
earlier flight control system designs with lower 
bandwidth requirements, it was satisfactory to use 
low-order analytic models which did not accurate- 
ly model the high-frequency rotor dynamics; with 
the recent introduction of high-performance, high- 
bandwidth control system specifications, it has be- 
come increasingly necessary to correctly model the 
coupled fuselage/rotor dynamic modes. It has long 
been known that flap dynamics introduce signifi- 
cant time delays into the rotor system, and more re- 
cently, Curtiss has shown that inclusion of the lag 
dynamics is important in the design of high perfor- 
mance control systems (Curtiss, 1986). Recent 
studies have explored the possibility of using rotor 
state feedback designs to damp blade motion (Ham, 
1983). An accurate understanding of the coupled 
fuselage/rotor dynamics is therefore important in 
rotorcraft control system design and stability analy- 
ses. 

Recent flight test experiments have shown 
that existing simulation models do not accurately 
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predict these high-frequency modes (B allin et al, 
1991, Kaplita et. al, 1989, and Kim et. al, 1990). 
These studies show significant differences between 
theory and experiment associated with the coupled 
rotor/body dynamics, especially in the frequency 
region dominated by the rota* lag motioa This re- 
search is therefore directed toward providing an im- 
proved understanding of the aeroeiastic and 
aeromechanical phenomena which determine the 
coupled rotor/body dynamics at hover. 

In order to gain physical insight into helicop- 
ter dynamics, development of linear models incor- 
porating coupled rotor/fuselage dynamics has long 
been a research objective. Past approaches to linear 
model development have included direct numerical 
perturbation of nonlinear simulations (D if tier, 
1988), identification of state-space stability and 
control matrix elements (Tischler, 1987), and ana- 
lytic derivation of linear equations of motion (Zhao 
and Curtiss, 1988). This study uniquely combines 
system identification methods with analytic model- 
ing techniques in order to investigate helicopter 
hover dynamics and to arrive at an improved linear 
model. The emphasis is on the high-frequency dy- 
namics of the coupled rotor/body motion. 

The identification study is carried out on 
flight test data from a Sikorsky H-53E helicopter at 
hover, using previously published data (Kaplita et 
al, 1987, and Mayo et. al, 1990). 

Research Objectives 

This paper describes an investigation into the 
response of articulated rotor helicopters to control 
inputs in hover. The goal is an improved under- 
standing of the coupled rotor/fuselage dynamics in 
hover directed toward a validated analytic Simula - 
tion model including high-frequency rotor/fuse- 
lage dynamics for use in stability analyses and 
high-performance control system design studies. 

Identification of linear, time-invariant state- 
space models representing high-order helicopter 
dynamics including main rotor degrees of freedom 
has long been an objective of engineers involved in 
rotorcraft simulation and control system design. 
The state and control matrix elements in an identi- 
fied state-space model can provide physical insight 


into system dynamics and can be used in combina- 
tion with mathematical modeling techniques to 
analyze differences between theory and experi- 
ment. 

State-space identification techniques have 
been applied to conventional fixed-wing aircraft 
with useful results. Since identification of state- 
space models using directly parametrized state and 
control matrix elements requires the estimation of 
a laige number of parameters, a reduced order mod- 
el is often used, assuming six degree-of-freedom 
rigid body dynamics and decoupling between the 
longitudinal and lateral axes. 

Identification of reduced order state-space 
models for rotorcraft have generally produced un- 
satisfactory results. The presence of the rota* pro- 
duces significant rotor/body coupling, requiring 
additional states to describe the high-frequency dy- 
namics, and also introduces significant interaxis 
coupling. The complete rotorcraft identification 
problem is therefore required to use a high- order, 
multi-input, multi-output model with as many as 
18 a more states. 

In order to avoid the inevitable problem of 
overparametriz ation which results when attempt- 
ing to identify a directly parametrized high-order 
helicopter model, this study uses an analytic model 
to generate state time histories. The model used in 
this study has been developed at Princeton using the 
Lagrangian formulation. It includes the coupled fu- 
sel age/rotor dynamics, main rota inflow, tail rotor 
thrust, and provides fa tail rotor inflow dynamics. 
It was analytically linearized about hover. This 
model provides a state-space description of the he- 
licopter at hover which is completely analytic and 
dependent only on an input set of physical parame- 
ters. A subset of these inputs are considered uncer- 
tain, and are to be estimated from flight test data. 
The flight-test derived parameter estimates can be 
used in combination with the mathematical fa- 
mulation to trace various physical aspects of 
coupled rota/body dynamics and thereby obtain 
physical insight The complete high-order model 
including rot or dynamics can be reasonably para- 
metrized by 15 or fewer physically meaningful in- 
put coefficients, resulting in a substantial reduction 
in the number of parameters to be estimated. 
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The framework of the identification ap- 
proach is the time-domain maximum likelihood 
methodology. The likelihood function is formu- 
lated assuming the presence of Gaussian measure- 
ment and process noise. The process noise may be 
nonwhite. The noise covariances as well as process 
noise dynamics may be parametrized. With Gaus- 
sian noise assumptions* the likelihood function be- 
comes the weighted least-square of the residual 
errors. The Kalman filter is the natural way to pro- 
duce these residuals for state- space dynamic sys- 
tems. 

The maximum likelihood estimate is ob- 
tained by fin ding the global maximum of the likeli- 
hood function. The parameters are nonlinearly 
related to the cost function and the resulting param- 
eter space is highly multimodal. Traditional func- 
tion optimization techniques based on gradient 
methods generally become trapped in local optima. 

The genetic algorithm is an alternative func- 
tion optimization approach which does not rely on 
the use of local gradient information. The genetic 
algorithm is an adaptive scheme* based on the anal- 
ogy with natural evolution* which efficiently 
searches a laige parameter space for the ’fittest’ 
solution to a given objective. This method has been 
demonstrated to be highly effective in obtaining the 
global maximum in a multimodal parameter space. 

The formulation of the system identification 
problem in the maximum likelihood framework 
leads to estimates of physical coefficients which 
have attractive statistical optimality properties and 
represent the best possible combination of physical 
coefficients necessary to match the given test data 
set. 

This identification methodology allows an 
assessment of model assumptions inherent in the 
mathematical model used to generate the state time 
histories. In this study, emphasis is placed on the 
frequency region associated with coupled rotor/fu- 
selage dynamics. In the frequency domain* the 
dominant feature in the rotor magnitude response is 
a notch characteristic produced by the presence of 
the in-plane blade degree of freedom. Using rotor 
blade constants derived through the identification 
procedure, rotor blade modeling assumptions may 
be examined, resulting in analytic model improve- 


ments. This study examines in detail the blade 
structural modeling assumption and investigates 
the effect erf accounting for blade flexibility effects 
generated by the presence of a large mechanical 
damper at the blade hinge. 

Analytic Model Description 

Research at Princeton has resulted in the de- 
velopment of a linearized rotor/body helicopter dy- 
namic model. The dynamic equations are 
formulated using a Lagrangian approach in order to 
capture all the important inertial coupling terms. 
The model includes rigid-body translation and 
rotation (pitch, roll, and yaw rates* longitudinal and 
lateral velocities)* rigid blade lag and flap multimo- 
dal coordinates, and main rotor cyclic dynamic in- 
flow. The controls are main rotor cyclic and pedals. 
The version of the model used in this study was ana- 
lytically linearized about the hover trim condition 
and does not include the collective degree of free- 
dom. 

Rotorcraft dynamics includes coupling be- 
tween the motion of the fuselage which is in rota- 
tional and translational motion relative to inertial 
space* and the motion of individual rotor blades. 
The final set of equations of motion are referenced 
to the body-fixed axis system which has its origin 
at the fuselage cento of gravity. In the Newtonian 
approach to modeling coupled rotor/fuselage equa- 
tions of motion, blade acceleration terms are first 
written referenced to the hub axis which is rotating 
at constant velocity; coordinate transformations are 
then used to obtain acceleration terms in the body- 
fixed frame. The complexity of the resulting accel- 
eration terms, combined with the number of 
degrees of freedom necessary to model rotor dy- 
namics properly, has led to the use of Lagrange’s 
equations for the derivation of the coupled rotor/ 
body model. 

The development of Lagrange’s equations 
proceeds from the evaluation of the Lagrangian, 
which requires only position and velocity terms in 
order to relate the system generalized forces to 
changes in the system kinetic and potential ener- 
gies. The generalized coordinates in Lagrange’s 
approach represent the degrees of freedom in the 
system and are chosen to correspond to the system 
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states. The kinetic energy term includes the motion 
of the fuselage and rotor blades, and the potential 
energy includes the gravitational potential energy 
of the fuselage and stored energy in the mechanical 
springs in the rotor system. Mechanical dampers 
are accounted for by use of the dissipation function. 
The generalized forces include aerodynamic forces 
due to fuselage and blade aerodynamics. Evalua- 
tion of the time and partial derivatives in the La- 
grangian can be time consuming for a high-order 
model and can be assigned to a symbolic manipula- 
tion program such as MACSYMA. 

Identification Methodology 

This paper describes an approach for identifi- 
cation of a coupled fuselage/rotor model for rotor- 
craft hover dynamics from flight test 
measurements. The identified model includes flap 
and lag degrees of freedom, main rotor inflow, and 
process and measurement noise disturbances. The 
process noise may be colored. The approach uses 
an analytically derived, linear time-invariant 
state-space model with literal coefficients which is 
parametrized in terms of aeromechanic al input co- 
efficients. The model order and structure may 
therefore be assumed to be determined by this ap- 
proach, and the system parameters are to be esti- 
mated from observations. The parameter 
estimation problem is formulated using the statisti- 
cal framework of maximum likelihood (ML) es- 
timation theory, thereby benefitting from known 
optimality properties of ML estimators. This dis- 
cussion first presents the parametrized dynamic 
model to be used in the identification methodology, 
and then describes the application of the maximum 
likelihood estimation approach to dynamic sys- 
tems. 

Model Parametrization 

The helicopter is modeled as a continuous- 
time dynamic system whose measurements are dis- 
cretely sampled as sensor outputs. Thus the 
identification algorithm is required to estimate con- 
tinuous-time model parameters from discrete sen- 
sor measurements. This continuous/discrete 
formulation is well known and is discussed by 


Ljung (1987). The linear time-invariant state 
equations are derived using the Lagrangian ap- 
proach, and are given by 

m = AJiem + BJtOMQ + F c (0Mt) ( 1 ) 

The model form accounts for the presence of pro- 
cess noise, where w(/) is assumed to be zero-mean 
white noise with unity spectral density. The contin- 
uous-time matrices, A c (0), 5/0), and F c (0), are 
parametrized by a vector of parameters, 0, which 
are to be estimated from observations. 

The observations are sampled at discrete 
time intervals, where 

Xm = COH m + G(0)vj(kT) 

t = kT, k = 0,1,2,... (2) 

and VjikT) are the disturbance effects at the 
sampled time intervals. 

For digital implementation of the identifica- 
tion algorithm, the continuous-time state equation 
given in Equation (1) is discretized using zero- cy- 
der hold. The input is assumed to be held constant 
over the sampling time interval, and the continu- 
ous-time state equation can then be integrated ana- 
lytically over the interval in order to obtain the 
discrete-time state equation. The zero-order hold 
discretization introduces a phase lag equivalent to 
one-half sample interval, which is taken into ac- 
count by advancing the control input by the corre- 
sponding one-half time interval. 

Eliminating time subscripts for simplicity, 
the discrete-time state-space equations are given 
by 

x(t+ 1) = MOMt) + B{Q)u{t) + F(0Mt) 
y(t) = CWMr) + GC0MO (J) 

This equation is now understood to be a discrete- 
time equation. Here, w(t) and v(r) are sequences 
of independent random variables with zero mean 
and unit covariance. 

Maximum Likelihood Formulation 

Let y^be a vector of observations which are 
supposed to be realizations of stochastic variables. 
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and let v( t ) be a multi-dimensional observation and y(t, 9) is generated using Equation (3) with the 
taken at time t: discrete-time Kalman filter formulation. 


Y N = [y(l).><2)....tfA0] 


The observations, Y N , depend on a vector of param- 
eters, 9, which are also considered to be random 
variables. The conditional probability density 
function for 9, given the observations. Y N , is then 
given by 


P(M*) - 


pwm : m 
P m 


(4) 


where p{9) is the prior distribution of the random 
parameter vector. A reasonable estimate for 0 can 
then be obtained by finding the value of 9 which 
maximizes the conditional density function given 
by Equation (4). With no prior knowledge of the 
distribution of 0, p(9) may be assumed to be uni- 
form. The best estimate for 9 is then obtained by 
maximizing the likelihood of obtaining the ob- 
servations. This leads to the ML, or maximum li- 
kelihood, estimator, given by 

9ml ■ argmax p(Y N \6) (5) 


For parametrized dynamical systems, with 
Gaussian noise assumptions, the maximum likeli- 
hood estimator has the form 


& ml ■ arg max p(Y N \0) 
o 


= arg max 

o 


t.o) - 

1 

jlog\A{Q)\ -^loglr 


(tf) 


where 

m — number of measurements 
e(t,0) « y(t) - %t.6) 

A(9) = Ee(G)e t (0) 


The Genetic Algorithm 

The evaluation of the likelihood function as 
presented in Equation (6) requires a search for the 
global maximum of the likelihood function over a 
multimodal parameter space whose contours are 
not known. Specifically, the identification method- 
ology has led to a function optimization problem 
where the performance measure is a highly nonlin- 
ear function of many parameters. The principal 
challenge facing the identification problem is the 
very large set of possible solutions and the presence 
of many local optima. Hill-climbing methods for 
function optimization based on finding local gradi- 
ents become trapped in local optima and are inade- 
quate far this problem. Genetic algorithms 
overcome these difficulties by efficiently searching 
the parameter space while preserving and incorpo- 
rating the best characteristics as the search prog- 
resses. 

The problem of function optimization can be 
addressed using the paradigm of adaptive systems, 
where some objective performance measure (the 
cost function) is to be maximized (i.e., adaptation 
occurs) in a partially known and perhaps changing 
environment The idea of artificial adaptive plans, 
based on an analogy with genetic evolution, was 
formally described by John Holland in the seventies 
and have recently become an important tool in 
function optimization and machine learning (Hol- 
land, 1975, and Goldberg, 1989). Holland’s artifi- 
cial adaptive plans have come to be known in recent 
literature as genetic algorithms. 

Genetic algorithms are based on ideas under- 
lying the process of evolution; i.e., natural selection 
and survival of the fittest. Using biological evolu- 
tion as an analogy, genetic algorithms maintain a 
population of candidate solutions, or ’individuals,’ 
whose characteristics evolve according to specific 
genetic operations in order to solve a given task in 
an optimal way. 

As a general overview, genetic algorithms 
have the following attributes which distinguish 
them from traditional hill-climbing optimization 
methods (Goldberg, 1989): 
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1 . GA’s work with a representation of the pa- 
rameter values rather than with the param- 
eters themselves. 

2. GA’s search from a population of points, 
not from a single point. 

3. GA’s use objective function information, 
not gradient information. 

4. GA’s use probabilistic transition rules, not 
deterministic ones. 

The genetic algorithm maintains a popula- 
tion of ’individuals’; i.e., possible solutions to the 
function optimization problem. In the context of 
the identification problem, each individual corre- 
sponds to a vector of parameters. The population 
of individuals ’evolves’ according to the rules of re- 
production and mutation analogous to those found 
in natural evolutionary processes, with the result 
that the population preserves those characteristics 
favoring the best solution to the cost function. 

The following steps were described by Hol- 
land (Holland, 1975) and contain the essentials 
properties of the the basic genetic algorithm 

1. Select one individual from the initial pop- 
ulation probabilistically, after assigning 
each individual a probability proportional 
to its observed performance. 

2. Copy the selected individual, then apply 
genetic operators to the copy to produce a 
new individual. 

3. Select a second individual from the popu- 
lation at random (all elements equally 
likely) and replace it by the new individual 
produced in step 2. 

4. Observe and record the performance of the 
new structure. 

5. Return to step 1. 

This deceptively simple set of instructions 
contains the ability to test large numbers of new 
combinations of individual characteristics and the 
ability to progressively exploit the best observed 
characteristics. It does so through the use of genetic 
operators. 


Genetic Operators 

Parent selection based on fitness, and the 
subsequent application of genetic operators to pro- 
duce new individuals are the steps by which the al- 
gorithm modifies the initial population and 
continually tests new combinations while main- 
taining those parameter sets which give high fit- 
ness. Each of these operations are performed 
probabilistically. 

The initial population of individuals is se- 
lected randomly with a uniform distribution over 
the defined parameter space. After one generation, 
parent individuals are selected randomly, with a 
probability which is proportional to the fitness as- 
signed to that individual. The selection procedure 
resembles spinning a roulette wheel whose circum- 
ference is divided into as many segments as there 
are individuals. The arc length of each segment is 
made proportional to the fitness value of the corre- 
sponding individual. Thus, the chance of choosing 
a given individual is uniformly random and yet pro- 
portional to its fitness. 

The genetic operations of crossover and 
mutation are then applied to the selected parent in- 
dividuals in order to introduce new characteristics 
into the population, enabling an efficient search for 
the optimal combination of parameters. 

The crossover operation involves a recom- 
bination of two selected individuals at a randomly 
selected point Thus the crossover operation pro- 
duces two new individuals, each of whom inheri t 
characteristics from both parents. 

The mutation operation involves a random 
alternation of an individual’s characteristic with a 
very low probability. This serves to introduce new 
information into the pool of structures and serves to 
guard against the possibility of becoming trapped in 
local optima. 

Genetic Coding 

Each individual is a candidate parameter set 
and is represented as a concatenation of individual 
parameters: 
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0 = [01 02 0.vl 


1. 


The genetic algorithm is illustrated in Figure 


In a digital implementation, each parameter 
Q l is encoded using a binary alphabet, and the indi- 
vidual is thus represented by a binary-valued 
string. The following specific coding scheme was 
suggested by Starer (Starer, 1990). 

Let each parameter 0, be bounded by Q imax 
and 0,. . If each parameter is coded in binary with 

a word length of /, then the interval | d imax , is 

discretized by 2 1 values. A representation of the pa- 
rameter 0 i can be obtained from the l-bit binary 
coding of 


mod 


- 0, )(2' - /) 

\ min J 


Q , - Q, 


To illustrate, let an individual represent a 
candidate parametrization where 


0 = [0i e 2 ) = [3 4.5] 


and bounds are given as 


0 10 0 
10 0 1 
1111 
0 10 0 
0 0 11 
10 11 


initial population 



parent selection 
based on fitness 


randomly selected 
crossover point 


crossover 


0 

1* 

1 

0 

1 

1 


0 

0 

1 random mutation 
0* 

0 

1 


J 

0 

0 

1 

0 

1 

1 


1 

0 

0 

1 new population 

1 
0 
1 


l<0 1 <4,2<0 2 <7,/ = 6 


Figure 1 The Genetic Algorithm 


The binary-valued string representing this candi- 
date vector is then 

Implicit Parallelism 

Genetic algorithms efficiently conduct a 
& binary - [101010 011111] search over a defined parameter space, converging 
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to a near-optimal solution. The basic unit of pro- 
cessed information in this genetic search is the 
schema, defined by Holland (1975). In the context 
of a digital implementation of genetic algorithms, 
a schema is a template specifying similarities at 
certain string positions. 

Thus, an individual is a string of binary dig- 
its, and the alphabet is composed of {0,7,#}, 
where # denotes ’don’t care’ (i.e., the value at this 
position has no effect on the performance measure- 
ment). As an example, an individual may be repre- 
sented as 


[0 0 1 1 1 0 1 1 0 0 0 1 0 ] 


A schema is a similarity template within this 
individual; so that this individual contains the sche- 
mata given by 

[0 0##101 100010 ] 


Given / positions, a single individual is an 
of 2 l distinct combinations, and an 
instance of 3 l distinct schemata. Further, a popula- 
tion of size N contains between 3 l and N3 l distinct 
schemata. Holland has shown that each schemata 
are evaluated and processed independently of the 
others, providing a tremendous computational le- 
verage on the number of function evaluations. 
Therefore, the use of genetic operators in the repro- 
ductive plan provides i) intrinsic parallelism in the 
testing and use of many schemata, and ii) compact 
storage and use of large amounts of information re- 
sulting from prior observations of schemata. 

The concept of implicit parallelism is funda- 
mental to the efficiency of genetic algorithms. 
Each schemata is processed and evaluated indepen- 
dently of other schema in the population; this pro- 
vides a tremendous computational leverage. A 
very weak lower bound states that for a population 
of (n) individuals, more than o(n 3 ) useful ’pieces’ 
of information is processed in each iteration (Gold- 
berg, 1989). 


An Example 

As an illustration of the genetic algorithm, 
consider the following example. 

/Uy) = 

3(1 - y ) 2 e -> 2 - < ~* +,)1 - 

/flg - / - x?'je~ y2 ~* 2 - L e -<*+» 2 -* 2 

The function surface is shown in Figure 2, 
along with the contour lines. This multimodal 
function has a global maximum at 
(7.557#, - 0.0093). 

A genetic algorithm was run on this function 
with a population size of 20. The initial guesses 
were chosen randomly, and were bounded as 
— 3 < x < 3, - 3 < y < 3. A binary code 
with wordlength of 8 was used, which means that 
both x and y were discretized by 256 points. An ex- 
haustive grid search under these conditions would 
involve evaluating 65536 possible points to find the 
global maximum 

Snapshots of the population distribution up to 
7 generations are shown in Figure 2. The snapshots 
show the population converging upon the global 
maximum; by the 7 th generation, most of the indi- 
viduals have converged on the maximum. The ge- 
netic algorithm in this case converges on 
( 1.5412 , — 0.0353) as the global optimum. 

This convergence has occurred after 7 gen- 
erations. With a population size of 20 individuals, 
this is 140 function evaluations as compared to the 
65536 necessary for the grid search. 

This relatively simple example serves to il- 
lustrate the ability of the genetic algorithm to find 
the optimum of a given function , using no gradient 
information , 

Analytic Model Validation 

The mathematical model is correlated with 
flight test data using nominal values for input coef- 
ficients. The correlation plots in Figure 3 show 
transfer function comparisons for pitch and roll 
axes. The data represent separate flights. In each 
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Sample Function 



Generation 1 


Generation 3 



Figure 2 Genetic Algorithm Example 
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case, the comparison is between the flight test rate 
gyro output and the model state. The comparison 
is made between 0.5 Hz (3.14 rad/sec) and 6 Hz 
(37.7 rad/sec) since the input signal was designed 
to cover this frequency range. The fuselage struc- 
tural bending modes are lightly damped and domi- 
nate the frequency above ~20 rad/sec. Therefore 
the identification procedure uses a bandpass filter 
with the upper cutoff frequency at 15.7 rad/sec. The 
frequency range of interest is therefore between 0.5 
Hz to 2.5 Hz (3.14 rad/sec to 15.7 rad/sec). 

The choice of physical coefficients used to 
parametrized dynamic model must allow adjust- 
ments to account for differences between test and 
theoretical responses using nominal physical input 
values. The gain differences at low frequencies, 
implying a mismatch in rigid body response, re- 
quires parametrization of the rigid body accelera- 
tion. The coupled fuselage/lagwise modes are a 
lightly damped pole-zero pair and create a notch- 
filter effect in the frequency response between 10 
-15 rad/sec. This frequency is near the -180 degree 
crossover, and a mismatch in this region adversely 
impacts the gain and phase margin calculations. 
Modeling the dynamics of this mode is important 
for control system design and stability analysis and 
will be the primary focus of modeling in this study. 


Validation Of Identification Procedure Using 
Simulated Data 

The maximum likelihood identification 
methodology for parametrized dynamic systems is 
validated first on a simulation with known parame- 
ters. These results demonstrate the feasibility of us- 
ing genetic algorithms to estimate physical 
coefficients from noisy data, and establish the pop- 
ulation size and crossover and mutation rates for 
this application. 

The simulation model is driven by flight test 
control inputs from the hovering condition. Main 
rotor pitch and roll cyclic and tail rotor pedals are 
all active, with primary excitation into roll cyclic. 
The output states used to form the cost function are 
pitch, roll, and yaw rates, and pitch and roll atti- 
tudes. No velocity information is necessary. 


Simulation Model Parametrization 

The model structure and parametrization was 
presented in Equations (1) through (3). The contin- 
uous-time state space model is analytically derived 
using the Lagrangian approach and using a vector 
of physical input coefficients, 0. For the purposes 
of this simulation study, the model structure has 
been augmented to include a first order time 
constant on process noise. The process noise dy- 
namics are to be parametrized and estimated from 
output data. 

The simulation model was parametrized as 
follows: 

aerodynamic coefficients: 
lift curve slope, a 

inflow equivalent cylinder height, hhnd 
inflow wake rigidity factor, wrf 

hover trim values: 

trim flap angle, fi 0 

trim main rotor pitch angle, ^ 

trim inflow velocity, v 0 

main rotor blade constants: 

lag damper constant, 
lag spring constant, 
flap spring constant, 

inertias: 

fuselage cross-moment, /** 
tail rotor: 

tail rotor thrust scale factor, K TR 
noise parameters: 

noise covariance ratio, NR 
process noise time constant, r 

Kalman filter theory allows optimal state estimates 
to be obtained in tire presence of state and measure- 
ment noise, where the Kalman gain is uniquely de- 
termined up to the ratio of process to measurement 
noise. The noise covariance estimate is therefore 
parametrized by the ratio of process to measure- 
ment noise. 

Genetic Algorithm Procedure 

The genetic algorithm was implemented us- 
ing a population size of 500 individuals; a crossover 


274 



««• 

Rad/Sec 


3 Model Validation 


F 5 



rate of 2/3; and a mutation rate of 1/1000. The pa- 
rameters were allowed to vary within 50 percent of 
the known simulation values. 



• s 4 ■ • • s u 


Iteration 

Figure 4 Best Likelihood Values 


The sensitivity of the cost function to the pa- 
rameter values vary widely. Therefore, as parame- 
ters begin to show convergence, the range of 
allowable values is progressively narrowed in order 
to demonstrate convergence for all parameters. 


The identification proceeds by r unning 
10-12 separate genetic algorithms simultaneously, 
where each algorithm begins with a new random 
number generator seed to select the initial guesses. 
Each set of runs therefore produces a scatter band 
of near optimal guesses for each parameter. The pa- 
rameters which influence the cost function most are 
identified most tightly. 

Figure 4 shows the progression of the best fit- 
ness values out of the population at each genera- 
tion. The results are shown in Figure 5. The solid 
line in each figure denotes the true value. 

The noise covariance ratio parameter cou- 
ples only very weakly to the cost function and dis- 
plays an almost random distribution until the 
physical coefficient estimates sufficiently con- 
verge. Therefore a two-step estimation procedure 
is required, where the noise ratio is allowed to re- 
main free until physical coefficients have con- 
verged. The physical coefficients are then fixed 
while the noise ratio is estimated. 

This methodology clearly demonstrates con- 
vergence. Twenty iterations of the genetic algo- 
rithm were run. Table 1 tabulates the parameter 
estimates. 


Table 1 Estimated Parameters, Simulation Study 


Parameters 

9o 

2s 

std 

lift curve slope, a 

5.73 

5.72 

3.98e-4 

inflow equivalent cylinder height, hhnd 

0.46 

0.46 

2.23e-4 

inflow wake rigidity factor, wrf 

2.0 

2.0 

1.34e-4 

trim flap angle, 

0.02 

0.02 

4.99e-7 

trim main rotor pitch angle, t 0 

0.05 

0.0497 

9.75e-6 

trim inflow velocity, v 0 

0.02 

0.0196 

2.61e-6 

lag damper constant, U g 

5.0 

4.978 

7.7e-3 

lag spring constant, 

75.0 

75.0 

7.06e-2 

flap spring constant, 

45.0 

44.92 

6.3e-3 

fuselage cross-moment of inertia, /** 

30,000 

30.035 

4.98 

tail rotor thrust factor, K TR 

1.0 

0.99 

9.35e-4 

covariance ratio, process/measurement, NR 

1.0 

0.97 

0.11 

process noise time constant, r 

-1.0 

-0.99 

1.8e-3 
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Flight Test Identification Results 

Data consistency checks ensure that errors in 
data collection do not interfere with the estimation 
procedure. The requirements for this step were 
minimal in this study, since this estimation method- 
ology requires only rate and attitude information. 
Consistency was checked by integrating accelera- 
tions and rates, and ensuring that sensor attitudes 
and rates match the integrated rates and attitudes. 

The flight test data was processed by 1) ap- 
plying a bandpass filter, and 2) decimating the data 
from 80 Hz to 8 Hz. The filter passband was from 
0 5 to 2.5 Hz (3. 1416 to 15.708 rad/sec). The lower 
bound corresponds to the beginning frequency of 
the frequency sweep input used to drive the system, 
and the upper bound is imposed to exclude the first 
fuselage bending mode at 3.4 Hz. 

The flight test identification parametrization 
was modified to reflect information available from 
comparison between test and theoretical responses 
generated from the analytic model using nominal 
parameter values. The parameter list used in flight 
test identification runs is shown in Table 2. The 
modifications are explained below. 

The parametrization of body inertias ac- 
counts for significant differences between theory 
and test in rigid body response, especially in the roll 
axis. Further, due to significant differences in 


cross-axis predictions, the roll and yaw rigid body 
responses could not be simultaneously satisfied. 
Therefore, yaw axis parameters were eliminated, 
and the identification scheme therefore attempts to 
fit pitch and roll responses only. This is permissible 
since for small motions about hover, yaw rate does 
not couple with main rotor cyclic multiblade coor- 
dinates and has no effect on pitch and roll responses 
in the rotor/body frequency region. 

The inflow equivalent cylinder height (hhnd) 
is related to the main rotor dynamic inflow time 
constant This parameter had no effect on the cost 
function in the bandpass frequency region used in 
this study. Therefore a quasistatic main rotor in- 
flow formulation was used and this parameter was 
dropped. 

The process noise dynamics, parametrized 
by a first order time constant, was also eliminated. 
This parameter is uniquely identifiable apart from 
the noise power ratio only if the time constant falls 
within the bandpass frequency range, and was 
found to have no effect cm the cost function. 

The identification run was carried out using 
flight test data from hover, with primary excitation 
into roll cyclic. The analytic model, parametrized 
as given in Table 2, was driven by main rotor pitch 
and roll cyclic and tail rotor pedal. The likelihood 
function was famed using pitch and roll rates only. 


Table 2 Estimated Parameters, Flight Test 


Parameters 

scale factor, fuselage roll moment of inertia, I x 
scale factor, fuselage pitch moment of inertia, I y 

lift curve slope, a 

inflow wake rigidity factor, wrf 

trim flap angle, fi 0 
trim main rotor pitch angle, ^ 
trim inflow velocity, v a 
lag damper constant, C g 
lag spring constant, K g 
flap spring constant, 
noise covariance ratio, NR 


A 


Qo 

std 

bounds 

nominal 

0.44 

0.011 

0.35-1.0 

1.0 

1.15 

0.033 

0.7- 1.3 

1.0 

8.4 

0.066 

5-10 

5.73 

8.0 

0.23 

2-11 

2.0 

0.162 

0.0013 

0.05-0.25 

0.0848 

0.0172 

0.00016 

0.005-0.15 

0.1304 

0.048 

0.0007 

0.01-0.1 

0.0613 

55 

0.10 

4-10 

9.5 

85.0 

0.735 

0-100 

0 

16 

1.34 

0-20 

0 

- 

- 

0.001-0.1 

_ 
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The initial choice of boundary limits on each 
parameter defines the parameter space to be 
searched in the identification algorithm. The 
bounds applied to each parameter are shown in 
Table 2; in each case, the bounds are chosen to in- 
clude the nominal value. 

Table 2 shows the identification results for 
flight test data. It was found that the noise ratio pa- 
rameter did not conveige while the remaining 
physical coefficients did, indicating that relative to 
the aeromechanical coefficients, noise powers af- 
fect the cost function only very weakly. 

The correlation with flight test data using the 
identified parameters is shown in Figure 6, where 
the roll axis response is correlated with the data set 
used in the identification, and the pitch axis re- 
sponse is an independent check. The roll axis cor- 
relation shows clear improvement in model 
correlation using identified coefficients. The low 
frequency gain prediction has been corrected 
through the inertia adjustment, and the notch in 
gain response due to the coupled lag/body response 
has been corrected. 

The differences between identified and nom- 
inal parameters can provide physical insight into 
rotor phenomena when analytic explanations can 
be found for parameter differences. The identified 
parameters for lift curve slope, a, and wake rigidity 
factor, wrf. \ have produced significant improvement 
in model response, indicating a possible require- 
ment for refinement of the aerodynamic theory 
used in the model. The identified parameters for 
main rotor spring and damping constants indicate 
necessary refinements in the prediction of frequen- 
cy and damping of blade motion. A model im- 
provement for blade in-plane dynamics is now 
presented. 

Modeling Blade Elasticity 

The identification procedure has resulted in 
estimated values for rotor blade spring and damp- 
ing parameters which are different from nominal 
values. The nominal mechanical damper value 
may be assumed to be known since it can be inde- 
pendently verified through available data. 


A procedure for modeling blade elasticity is 
presented which accurately accounts for differ- 
ences between nominal and estimated values for in- 
plane motion frequency and damping. The method 
of assumed modes is used to model the case of a 
flexible beam with damper and spring constraints. 
This procedure is first demonstrated on a nonrotat- 
ing beam, for which an exact solution can be ob- 
tained. The method of assumed modes will be 
shown to be a good approximation of the exact solu- 
tion. This approximate solution can then be used in 
the flexible beam analysis in the analytic hover he- 
licopter model. The beam formulations for both ro- 
tating and nonrotating blades with both spring and 
damper constraints at the root is given in detail in 
Appendices A and B. 

Approximate solution methods such as the 
method of assumed modes display convergence to- 
ward the analytic solution as more assumed mode 
shapes are added to the set of basis functions. The 
first approach to the lagwise bending problem was 
to use increasing numbers of mode shapes that ful- 
filled the boundary conditions for a hinged beam. 
However, with this approach, convergence was not 
achieved after even after using 5 assumed modes. 
In order to avoid using an unacceptably large num- 
ber of basis polynomials in the model, an alterna- 
tive approach using a combination of modes that 
satisfy hinged and cantilever boundary conditions 
was used. 

Figure 7 illustrates the assumed modes solu- 
tion method using both the nonrotating and rotating 
beam formulations. For a nonrotating beam with 
spring and damper constraints, an exact expression 
for the beam eigenvalues is available and is given 
in detail in Appendix B. The analytic eigenvalue 
equation is solved numerically. In this case, the 
root finding problem was converted into a function 
optimization problem and solved using the genetic 
algorithm. This solution to the exact formulation is 
shown against approximate solutions in Figure 7. 
The approximate solution using the Lagrangian ap- 
proach, when using only basis functions which ful- 
fill hinged beam boundary conditions, approach the 
exact solution slowly. With 4 hinged basis polyno- 
mials, the solution has not yet converged. Howev- 
er, the assumed modes approach with only one 
hinged plus one cantilever mode shapes matches 
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the analytic solution exactly. Convergence is dem- 
onstrated by the fact that addition of either hinged 
or cantilever mode shapes do not further change the 
eigenvalue solution. 

Figure 7 then shows the convergence of the 
approximate solution for the rotating beam, for 
which there exists no known exact solution. Here, 
the sum erf 2 hinged plus 2 cantilever modes is near 
convergence. The addition of either one more 
hinged or one more cantilever mode does not 
change the solution appreciably. The combination 
of 2 hinged plus 2 cantilever modes is chosen for 
model development as a good compromise between 
model order and accuracy of solution. 


iw 

Q 


o_ 

Q 

Figure 8 Rotating Frame Lag Roots 


Rigid , Identified - 
Elastic 




Rigid, Nominal 


Figure 8 shows the location of the rotating 
frame lag mode eigenvalues. The elastic blade 
model using two hinged and two cantilever mode 
shapes is used to show the progression of the root 
location as damper value is increased from zero to 
the nominal value. The predicted root location for 
the elastic model with the nominal damper constant 
agrees reasonably well with the predicted location 
for the rigid blade model using a fictitious spring 
and using identified spring and damper constants. 
The rigid blade model using nominal damper 
constant only (no spring) predicts a much higher 
damping and lower frequency than is indicated by 
test data. 


Conclusions 

An analytically derived linear model of 
coupled rotor/body dynamics at hover has been val- 
idated against flight test data. 

The analytic model with literal coefficients 
has been parametrized using 11 physically mean- 
ingful coefficients, including noise covariances. 
This model has been used to formulate a multi-in- 
put, multi-output likelihood function in the time 
domain. The analytic model is used to generate the 
state time histories. Only body rates are necessary 
in the cost function. 

The likelihood function is globally maxi- 
mized using the genetic algorithm approach, result- 
ing in statistically optimal maximum likelihood 
parameter estimates. 

The estimated parameters indicate that lag 
mode damping in flight is approximately one-half 
of the value expected from rigid blades. 

The correct analytic prediction for lagwise 
motion is obtained using an elastic blade formula- 
tion. The flexible blade model was formulated us- 
ing a normal mode approach and checked using the 
closed form solution for a nonrotating beam. The 
convergence results using assumed mode shapes in- 
dicate that the correct lagwise bending mode 
shapes are obtained using a combination of cantile- 
ver and hinged assumed modes. 
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Appendix A. Modeling Blade Elasticity 

Equation (A.l) gives the in-plane bending 
equation for a rotating beam. The derivation can be 
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found in Bramwell (1976), and in Johnson (1980). 
This partial differential equation relates the mo- 
ments due to the inertial, centrifugal, and aerody- 
namic forces to the moment expression from 
engineering beam theory. 



All quantities are understood to refer to lagwise 
bending motion. Here, G(r) is the centrifugal ten- 
sion force at a point at a distance r from the hub cen- 
ter, £ is the modulus of elasticity, / is the lagwise 
area moment, and Q is the rotor rotational velocity. 


I x i ei X q " r2 *x g X q ” dx < * >n + 


X& - Q2q " = ° 


(A.3) 


Multiply Equation (A.3) by <p m and integrate from 

% < x < R, or ?< x < 1 where F is understood 
l\ 

to be a nondimensional offset value. 

This gives 






+ R 4 Yjq* ~ "*t>*4><4x = 0 {A. 4) 


The boundary conditions for a hinged blade 
are: 


Integrating each term by parts, the first term gives 


At the hinge: 

Y{e) = 0 

£/4-t = moment = 0 

dr 2 


At the tip: 


£/ tt = 0 

dr 2 

= shear force = 0 


There is no known analytic solution for 
Equation (A.l) due to the presence of the centrifu- 
gal term. A solution based on the method of as- 
sumed modes is presented. 


Let the lagwise displacement be of the form 


1 r V 

" I L " Jr 

r i 1 1 

- + | X 4rM<<Pn4>n4x 

= X <t>'^t> m q n + I X q£ l( S>H<Prr£x (A_5) 
n J n 


Equation (A.4) was obtained using the boundary 
conditions for the hinged blade, along with the end 
constraint imposed by the damper, which is given 
by 


Y(x.t) - R^<P*(x)q„(t) (A2) 

n 


where R = blade length. This solution method fol- 
lows the method of separation of variables. <p n (x) 
are a sequence of functions, not necessarily ortho- 
gonal, which approximate the expected blade shape 
and which satisfy the blade boundary conditions. 

Substituting into Equation (A.l), 


E! m - 
El *■!„. 


- D 


d 2 Y 


-DR 

dx dt I 


where D = damping constant. 
Similarly, the second term gives 



tPm^Gip'ndx » 
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1 

-K 2 X«»j 


(A. 6) 


<p{x) = x 6 - 2X 5 - + ^x 3 + x 


Using Equations (A.4) through (A.6), 


cantilever-free: 



q£l<p' n <pjx 4- RD ^ 

/i 


+ i? 2 



Gtpntfrrrtdx 


+ i? 4 ^(<7„ - £ 2 <? n ) nup^dx = 0 


- x 4 - 4r* + dr 2 
<pM » ^ - Kk 3 + 2Qr 

Since these polynomials meet boundary 
conditions at x=0 and at x~L and the blade for- 
mulation is integrated from x = I to x=l , the basis 
polynomials are transformed to new coordinates, 
where 

x = (1 - ?)x + 


To evaluate this, nandimensionalize by m£> 2 R 4 and 
collect terms, which results in 

AnmQ n 4 * ^nmQ a ^rw^Jn ^ 

where 


With this coordinate transformation, the new set of 
polynomials, which now fulfill the necessary 
boundary conditions at the hinge offset and at the 
blade tip, are now 

hinged-free: 


i 


__ RD£2 j ' j ' 

Dnm 

l 


Bnm ~ j \m&R* < t >n< t >m + *** 


Basis Functions For Assumed Mode Shapes 


<t>(x) = X - ? 

4>{x) = 1.4& 6 - 3.33X 5 - 0.12* 4 + Alx 3 - 
0.79X 2 + 1.12* - 0.07 


cantilever-free: 

<t>(x) = l.3x* - 5.2x* + 7.S* 2 - 0.92* + 0.03 
<p(x) = 1.39* 5 - 0.44* 4 - 12.1 lx 3 + 

25.10* 2 - 3.03x + 0.09 


Polynomials are used as the basis functions, 
<p„(x). Two sets of polynomials, meeting the neces- 
sary boundary conditions for hinged-free and can- 
tilever-free beams, were used in this study. They 
are: 

hinged-free: 


Appendix B. Exact Equations Of Motion For 
A Nonrotating Beam 

The modal analysis assumes that the beam 
displacement is written as a sum of modal displace- 
ments: 


/v KCx.f) = RY<P„(x)q,W 

<p(x) « X , 
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To find the exact analytic solution in the case 
of root constraint with both spring and damper, note 
that the boundary conditions are given by 

<p{0) = 0 

<t>'X0) = [§ + &$] 4>X0) 

4>"U) = 0 
<P"‘U) = 0 

where K and D are spring and damper constants and 
all quantities are understood to refer to lagwise mo- 
tion and are defined as in Appendix A. 

These boundary conditions are satisfied by 
writing the mode summation equation as 

<p{x) = fcW + [§ + a 

where <pf(x) and <p c (x) refer to hinged and cantile- 
ver mode shapes. 

The hinged end mode shape solutions are 
given by 

(pfix) = cos(A ) sinh(Ax) + cosh(A) sin(Ax) 

(ppiO ) = 0 

<Pf(0) = A [cos(A) + cosh(A )] 
fpi.0) = 0 

<Pf(1) = A 2 [cosiA) sinh(A) - coshiA) sin(A)] 

<p F (l) = A 3 [cosiA) cosh(A) - coshiA) cosiA)]] 

The cantilever mode shape solutions are giv- 
en by 

<PcW = ( sin(A ) - sinh(A))(sin(Ax) - sinh(Ax)) + 

(cos(A) + cosh(A))(cos(Ar) - cosh(Ar)) 

<Pc(0) = 0 
<PcW) = 0 

<£c(0) = - 1A 2 [cos(A) + cosh(A) ] 

<t>c£l) = - A 2 [1 + cosh(A) cos(A)] 

<f>c(l) = A J [t««(A) - sinh(A))( - cos(A) - cosKA)) H 
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(casfA) + cosh(A))(sin(A) - $wA(A)) ] 
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